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1. Introduction 

The impressive success of the miniaturization of technical devices down to 
submicrometric sizes has led to numerous practically applicable devices which are 
used for the manipulation and transport fluids. Such microfiuidic devices include 
micromixers, DNA amplifiers, chromatography systems, or chemical reactors. However, 
in order to fully understand the properties of such devices and to optimise them for 
efficient usage, some fundamental questions in physics including the behavior of single 
molecules in fluid flow or the validity of the no-slip boundary condition have to be 
answered [ll, y] • While for gas flows at large Knudsen number rarefaction effects that 
can result in a slippage of the flow at surfaces are expectecd, for liquids one would 
naively assume that its velocity close to a surface always corresponds to the actual 
velocity of the surface itself. This assumption is called the no-slip boundary condition 
and can be accounted for as one of the generally accepted fundamental concepts of fluid 
mechanics. However, this concept was not always well accepted. Some centuries ago 
there were long debates about the velocity of a Newtonian liquid close to a surface 
and the acceptance of the no-slip boundary condition was mostly due to the fact that 
no experimental violations could be found, i.e., a so-called boundary slip could not be 
detected. 

Recently, well controlled experiments have shown a violation of the no-slip boundar 
condition in sub-micron sized geometries. Since then, mostly experimental [2|-|9| 



but also theoretical works ^0|, lUj, as well as computer simulations |12l-ll6| have been 
performed to improve our understanding of slippage. An (apparent) violation of the 
no-slip boundary condition can be explained by the complex behavior of a fluid close to 
a solid interface which involves the interplay of many physical and chemical properties. 
These include the wettability of the solid, the shear rate or flow velocity, the bulk 
pressure, the surface charge, the surface roughness, as well as impurities, dissolved 
gas, or bubbles attached to surfaces. Due to the large number of different parameters, a 
signiflcant dispersion of the results can be observed for almost similar systems [2, 113]. For 



example, observed slip lengths vary between a few nanometres [18| and micrometres ^^ 
and while some authors flnd a dependence of the slip on the flow velocity [g, [l9|, yj , others 
do not p, |j] . Extensive reviews of the slip phenomenon have recently been published 
by Lauga et al. J2|, Neto et al. [IT], as well as Bocquet and Barrat 20 . 

A boundary slip is typically quantifled by the so-called slip length 6 - a concept that 
was already proposed by Navier in 1823. He introduced a boundary condition where 
the fluid velocity at a surface is proportional to the shear rate at the surface |2ll] (at 
X = xo), i.e., 

«.(.») ^^^. (1) 

In other words, the slip length b can be deflned as the distance from the surface where 
the flow velocity vanishes. 

During the last few years, the substantial scientiflc research invested in the slip 
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phenomenon has lead to a more clear picture which can be summarized as follows: 
one can argue that many surprising results published were only due to artefacts or 
misinterpretation of experiments. In general, there seems to be an agreement within 
the community that the no-slip boundary condition is a valid assumption for smooth 



hydrophilic surfaces down to contact [9|, l22|, |23|. Slip lengths larger than a few 



nanometres can usually be referred to as "apparent slip" and are often caused by 
experimental artefacts or not fully understood effects. Examples for those can be surface 
structures, variations in the local contact angle at the surface or dissolved gases forming 
nano- or microscale bubbles on the surface. It is also possible to utilize such surface 
properties in order to generate bubbles in a controlled way so as to add slippery surfaces 
to the flow channel and thereby increase the slip length and reduce the viscous friction. 
An example that has lately gathered a lot of interest are surfaces that are patterned with 



arrays of holes where gas bubbles are attached to the holes [2J-|29|. Although reasonably 
flat bubbles lead to increased effective slip lengths, it has recently been observed that 
bubbles strongly protruding to the flow channel may actually have an opposite effect, 
i.e., they may reduce the slip and increase the hydrodynamic dra g \ 2m . Numerical 



studies have predicted even a possibility of negative slip lengths [24, l25|, |28[. The 
enhanced friction can be understood with help of the additional roughness produced 
by the bubbles which can overcome the lubricating effect of the stress-free liquid-gas 
interfaces. The numerical results of negative slip lengths were recently verified in 



theoretical investigation of Davis and Lauga [26 . 

In all studies concerning slip fiow over a bubble mattress we are aware of, the 
mattress under consideration has been built up of evenly sized and equally shaped 
bubbles. However, as real configurations inevitably include some kind of imperfections, 
it would also be interesting to know whether the disorder in the mattress affects the 
detected slip behavior. Therefore, in this paper a case where the mattress is made up of 
bubbles of two different sizes is considered. In practice, such a situation is achieved by 
structuring the surface with holes that have different radii. This structuring leads, for 
a given bulk pressure, to a bubble configuration where bubbles located in larger holes 
protrude stronger to the flow channel than those in the smaller ones. In particular, we 
show that the bubble-size distribution can significantly affect the slip behavior and, e.g., 
decrease the negative slip caused by strongly protruding bubbles. While this reduction 
may hamper experimental observations of negative slip, it may on the other hand be 
utilized in designing bubble surfaces aiming to maximal flow slippage. 

2. Lattice Boltzmann simulations to investigate boundary slip 

The simulation method used to study microfluidic devices has to be chosen carefully. 
Molecular dynamics simulations (MD) are the best choice to simulate the fluid-wall 
interaction, but the presently available computer power is not sufficient to simulate 
length and time scales necessary to achieve orders of magnitude which are relevant 
for experiments. However, boundary slip with a slip length b of the order of many 
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molecular diameters has been studied with molecular dynamics simulations by various 



authors 0, Q,!!!, S, 



31 



The current contribution focuses on numerical investigations of the slip phenomenon 
by means of lattice Boltzmann simulations. It should be noticed that while a large 
number of groups utilizes the lattice Boltzmann technique to investigate microfiuidic 
problems, only a very small number of researchers is actually applying the method to 
study shppage. This is surprising since mesoscopic simulation methods offer a closer 
relation to experimentally relevant time and length scales than microscopic techniques 
such as molecular dynamics, even though the interactions between fluids and surfaces 
have to be described on a mesoscopic scale in the lattice Boltzmann method. 

Within the lattice Boltzmann method a popular approach is to introduce slip by 
generalizing the no-slip bounce-back boundary conditions in order to allow specular 



reflections with a given probabihty [12|, |l6|, l32|, |33| , or to apply diffuse scattering 
36|. It has been shown by Guo et al. that these approaches are virtually equivalent 
Another possibility is to modify the fluid viscosity due to local density variations in 
order to model slip [38]. In both cases, the parameters determining the properties at 
the boundaries are "artificial" parameters and they do not have any obvious physical 
meaning. Therefore, these parameters are not easily mappable to experimentally 
available quantities. We model the interaction between hydrophobic channel walls and 
the fluid by means of a multiphase lattice Boltzmann model. Our approach overcomes 
this problem by applying a mesoscopic force between the walls and the fluid. A similar 
approach is used by Zhu et al. 39|], Benzi et al. 40|, and Zhang et al. 4l|]- This force 



applied at the boundary can be linked to the contact angle which is commonly used by 
experimentalists to quantitatively describe the wettability of materials 42l-l45| . While 
we are not aware of further lattice Boltzmann simulations to study the flow over a 
bubble mattress, a number of authors has applied various lattice Boltzmann multiphase 
and multicomponent models to study the properties of droplets on chemically patterned 
and superhydrophobic surfaces [461450 . 



The lattice Boltzmann method is based on the Boltzmann kinetic equation 



d 

^ + ^-^^ 



/(x,u,t) = n, 



(2) 



which is discretized on a lattice. The Boltzmann equation describes the evolution of 
the single particle probability density /(x, u, t), where x is the position, u the velocity, 
and t the time. The derivatives on the left-hand side represent propagation of particles 
in phase space whereas the collision operator Q, takes into account molecular collisions. 
When discretized, to represent the correct hydrodynamics, the collision operator should 
conserve mass and momentum, and it should ensure sufficient isotropy and be Galilean 
invariant. By performing a Chapman- Enskog analysis, it can be shown that such a 



collision operator Q leads to flow behavior following the Navier-Stokes equation 51 



In the lattice Boltzmann method the time t, the position x, and the velocity u are 
discretized leading to a discretized version of Eq. [2) 



/i(x + Cj, t + 1) - /i(x, t) = ni, 



0,1, 



B. 



(3) 
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Above, /j(x, t) indicates the amount of fluid at site x at time step t with velocity Cj,. 
Our simulations are performed on a three dimensional lattice with B = 19 discrete 
velocities (the so-called D3Q19 model). For the collision operator Qi we choose the 
Bhatnagar-Gross-Krook (BGK) form 52 



1 



a = — (/.(x,t)-/r(u(x,t),f?(x,t))) , 



r 



where r is the mean collision time that determines the kinematic viscosity 
2r- 1 



u 



6 



(4) 



(5) 



of the fluid. In this study the relaxation time r is kept constant at value 1.0. Due 
to the collisions, the system relaxes towards an equilibrium distribution f^'^ which can 
be derived imposing restrictions on the microscopic processes, such as explicit mass 
and momentum conservation. In our implementation we choose for the equilibrium 
distribution function 



fr = oq 



"l , Ci-U ^ (Ci-U)2 


u^ 


r oi ' 2cf 


'^cl 



(6) 



which is a polynomial expansion of the Maxwell distribution where Q are the lattice 
weights resulting from the velocity space discretization, and Cs = l/-\/3 is the speed of 
sound for the D3Q19 lattice. The hydrodynamic quantities are obtained as moments of 
the single-particle distribution function /j(x, t). For example, the density at lattice site 
X is 



i 

and the macroscopic velocity u(x, t) is obtained from 

£»(x,t)u(x,t) = ^/i(x,t)Ci. 



(7) 



(8) 



Mean-field interactions between fluid particles are introduced by following the work 
of Shan and Chen, as a mean- field body force between nearest neighbours 53|, l54| . 



G'6^(r)^CiV^(r + Ci)ci 



(9) 



where ijj = 1 — exp {—q/qq) is an effective mass, G^ tunes the strength of the interaction, 
and ^0 is a reference density. This force term leads to a non-ideal equation of state 
with pressure P = cj.g + ^c^Gbip'^, and it enables simulations of liquid-vapor systems 
with surface tension. To model the wetting behavior at fluid-solid surfaces, a similar 
interaction is added between the fluid and solid phases, and the contact angle is 



tuned by setting a density value g^, at the boundaries [55|, |45|. Additionally, we 
apply mid-grid bounce back boundary conditions between the fluid and the surface 
which assures vanishing velocities at solid surfaces. Here, a distribution function that 
would be advected into a solid node is simply reversed and advected into the opposite 



direction 51 



From molecular dynamics simulations it is known that the fluid-wall interactions 
causing a slip phenomenon usually take place within a few molecular layers of the liquid 
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along the boundary surface [l5|, |lj, l30|, lZ|. Our coarse-grained fluid wall interaction acts 
on the length scale of one lattice constant and does not take the molecular details into 
account. Therefore, coarse-grained implementations based on the lattice Boltzmann 
method are only able to reproduce an averaged effect of the interaction and cannot fully 
resolve the correct flow profile very close to the wall and below the resolution of a single 
lattice spacing. However, in order to understand the infiuence of the hydrophobicity on 
experimentally observed apparent slip, it is sufficient to investigate the fiow behavior on 
more macroscopic scales as they are accessible for experimental investigation. Coarse- 
grained interaction models could be improved by a direct mapping of data obtained 
from MD simulations to the coupling constant Gb allowing a direct comparison of the 

Similar approaches are 



45 



influence of liquid-wall interactions on the detected slip 

known from quantitative comparisons of lattice Boltzmann and molecular dynamics 

simulations in the literature [56|, l57 . 

In recent years we studied the influence of a number of parameters on an apparent 
boundary slip using the method described above. Here, we shortly review the main 
conclusions of our previous flndings. A more comprehensive review can be found in 58 . 



In 45 



we showed that our mesoscopic approach is able to reach the small 
flow velocities of known experiments and reproduces results from experiments and 
other computer simulations, namely an increase of the slip with increasing liquid-solid 
interactions, the slip being independent of the flow velocity, and a decreasing slip with 
increasing bulk pressure 45 . 

If typical length scales of the system are comparable to the scale of surface 
roughness, the effect of roughness cannot be neglected anymore. The influence of surface 
variations on the slip length b has been investigated by numerous authors. It was demon- 
strated by Richardson that roughness leads to higher drag forces and thus to no-slip on 
macroscopic scales 59|. An experimental conflrmation was later presented by McHale 



and Newton [60|]. Sbragaglia et al. applied the LB method to simulate fluids in the 



vicinity of microstructured hydrophobic surfaces [55|, Al-Zoubi et al. demonstrated 



that the LB method is well applicable to reproduce known flow patterns in sinusoidal 



channels 



61 



and Varnik et al. 



62 



63 



have shown that even in small geometries rough 
channel surfaces can cause flow to become turbulent. Recently, we presented the idea 
of an effective wall position at which the no-slip boundary condition holds for rough 
channel surfaces 6J]. We investigated the influence of different types of roughness on 
the position of the effective boundary hes- Further, we have shown how the effective 
boundary depends on the distribution of the roughness elements and how roughness and 
hydrophobicity interact with each other 65|- We were also able to simulate flow over 



surfaces generated from AFM data of gold coated glass used in microflow experiments 



by Vinogradova and Yakubov 66|]. We found that the height distribution of such 
a surface is approximately Gaussian and that a randomly arranged surface with a 
similar distribution gives the same result for the position of the effective boundary 
although in this case the heights are not correlated 6^. Currently, we model AFM 
based measurements to determine lubrication forces in order to probe the fundamental 
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concept of boundary conditions leading to a slippage or a shift of the effective boundary 



position |67, 



68 



A natural continuation of our previous works on roughness induced apparent 
boundary slip is the analysis of flow along superhydrophobic surfaces as presented in 
this article. It has been recently predicted [69| and experimentally reported 70|] that 
the so-called Fakir effect or Cassie state considerably amplifies boundary shppage on 
superhydrophobic surfaces. Using highly rough hydrophobic surfaces such a situation 
can be achieved. Instead of entering the area between the rough surface elements, the 
liquid remains at the top of the roughness and traps air in the interstices. Thus, a 
very small liquid-solid contact area is generated. In this article we quantify the slip in 
such systems using Couette and Poiseuille flow, where the flow is confined between two 
parallel walls separated by a distance 2d. In our simulations, the lower wall is static 
and is patterned with holes or grooves to which vapor bubbles are trapped. The upper 
one is smooth and can be driven with shear velocity Uq in the case of a Couette setup. 
The system boundaries are periodic and a unit cell of the regular array is included in 
a simulation. To trap bubbles to holes, some heterogeneity is needed at the edges of 
holes in order to pin the contact line. To this end, we use different wettabilities for 
boundaries in contact with the main channel and with the hole. For a Couette setup, 
the effective slip length can be calculated from the shear stress a = fidu/dz acting on 
the upper wall, which is obtained from the no-slip boundary condition imposed at the 
fluid-solid boundaries. Thus, the effective slip length reads as fe = fiuo/a — 2d, where n 
is the dynamic viscosity of the liquid. In case of a Poiseuille flow, a driving body force 
is applied effectively generating a pressure gradient ^. Assuming Navier's boundary 
condition, b is measured by fitting the theoretical velocity profile as given by 



uAx 
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26 



(10) 



2(^2 ■ d + b ' d + b_ 

to the simulated data via the slip length b. This flow profile assumes the lower boundary 
to have a slip while on the upper wall a no-slip boundary condition is applied. 




Figure 1. Definition of the protrusion angle tp. 



The parameter Gb is chosen such that the density ratio between liquid and gas is 
between 22 and 44. This ratio is too small for a realistic description of gas bubbles in a 
liquid. Also, the interface between both phases is of finite width causing the resistance 
in the vapor phase. However, these limitations of multiphase LB models do not influence 
the qualitative insight obtained from the simulations. 
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In equivalence to the definition of the contact angle, which is not suitable for bubble 
matresses as studied in this paper, we utilize the so-called protrusion angle cp to classify 
the bubbles in the simulated systems. The protrusion angle is defined by measuring 
the angle between the surface and a tangent at the bubble surface (see Fig. [1]). The 
protrusion angle can be varied in the simulations by changing the bulk pressure of the 
liquid phase. 

3. Results 

3.1. Surfaces patterned with microgrooves 

A surface structure that has been studied both experimentally and theoretically in 
great detail consists of longitudinal grooves etched in a smooth surface. If such grooves 
are hydrophobic, they can be filled with vapor forming cylindrical bubbles that are 
surrounded by liquid. The surface between individual grooves is hydrophilic causing 
an effective pinning of the three-phase contact line. In simulations, such hydrophilic 
surfaces can be obtained by setting a virtual surface density of the same order as the 
liquid density when computing the fluid-surface interactions as given by Eq. |9l The 
hydrophobic interactions are generated in a similar fashion by choosing a substantially 
smaller value for the virtual surface density. Fig. |2] shows a 3D visualization of the 
simulated unit cell. The length of the system is L = 40 lattice units and the width of 
the groove is defined by c resulting in a distance between two grooves given by {L — c) 
lattice units. The height of the channel is d = 135 lattice units in x-direction, where 
the top of the system is defined by a smooth no-slip surface. The flow is generated by 
a constant body force acting in 2;-direction effectively resulting in a Poiseuille profile. 
Apart from the top and bottom plane, all boundaries are periodic. 




Figure 2. 3D visualization of a grooved surface. The width of the unit ceh is L and the width 
of a groove is defined by c. 

I — I 

A similar system was investigated in the experiments by Tsai et al. [29;]. They 
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studied the flow field over gas filled grooves oriented in flow direction by using micro 
particle image velocimetry (micro PIV). In 29| they presented measured flow profiles and 
the corresponding effective slip lengths. The behavior of the slip length in dependence 
on different geometrical parameters has been investigated and it was concluded that 
the obtained values for slip length cannot be explained by the theory of Philip 71 
which provides an analytical model for slip flow on striped full-slip / no-slip surfaces. 
Their explanation was that due to the meniscus forming the bubble surface the no- 
shear (i.e. full-slip) parts develop an additional drag, similar to the effect reported by 
Richardson who showed that a rough no-shear surface creates sufficient drag to obtain 



macroscopically a no-slip boundary 59 
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Figure 3. Flow profile perpendicular to the grooves for flow parallel (left) and perpendicular 
(right) to the grooves at different x positions over the surface. The groove is located in the 
centre of the system between 0.25L and 0.75L. In the left panel a strong increase of the velocity 
in close vicinity to the bubble can be seen that relaxes within 0.1 6d to an almost undisturbed 
flow. The right plot demonstrates that flow perpendicular to the grooves causes more significant 
disturbances, which get still damped rather quickly. 



Following the work of Tsai et al. we first present simulations of flow over grooves 
oriented in flow direction. The width of the groove is chosen to be c/L = 0.5. The 
protrusion angle in the experiments and in the simulations is chosen to be y9 = —10°, 
i.e., the liquid-gas interface forms a dimple rather than a bubble. The meniscus increases 
the vapor covered area (no-shear surface) but also introduces roughness that leads to 
a higher friction and therefore a reduced slip. In Fig. [3] we show the flow profile over 
one groove for different distances from the surface. At a very close distance x = O.OAd 
the bubble has a strong influence on the flow, since the velocity in the centre of the 
groove is more than twice as large as the velocity in the centre of the void space. This 
effect is damped very quickly, i.e., at a; = O.lQd the velocity increase in the area above 
the bubble is less than 10% compared to the undisturbed case. The strong damping 
allows the treatment of such a grooved surface as a surface with an effective boundary 
position. The flow profiles are consistent with the results of Tsai et al. who also found 
a strong increase of the velocity above the grooves that is damped strongly further in 
the bulk. Further, the measured and the simulated slip lengths b/L are of the order of 
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half the theoretical prediction by Philip 7l| . The reason for this is the assumption of a 



flat surface with stripes of no-slip and no-shear, while in the case of vapor filled grooves 
a meniscus is formed introducing an additional roughness that reduces the slip of the 
fluid over the surface. 

In order to compare our results to the theoretical results of Lauga and Davis [26i], we 
also consider grooves which are oriented perpendicular to the flow. In the right panel of 
Fig. ini we show the flow profile in the vicinity of a perpendicular groove as it is assumed 
by the theory. Interestingly, here the flow velocity does not increase as strongly as in 
the aligned case above the groove but right after it. Again, a deceleration in front of the 
groove can be observed and the flow profiles are not symmetric towards the centre of a 
unit cell which is caused by the driving direction of the flow. The disturbances caused 
by the bubble are not damped as quickly as in the case of longitudinal grooves. For 
example, even at x = O.lGd the flow velocity changes between 0.45fmax and 0.59fmax- 
Further the disturbance close to the boundary is much stronger. At x = O.OAd one 
can see a strong decrease of the velocity at the beginning of the groove and a strong 
increase at the end. The reason is the dimple shape of the meniscus. This leads to 
a compression of the streamlines at the end of the meniscus which is equivalent to an 
increased velocity. The dependence of the slip length on the orientation is consistent 
with the work of Bazant et al. who have shown for fiat stripes that a generalized tensor 



form would be required to describe the surface properly [72 . 

The shape of the meniscus has a significant influence on the slip length. The above 
mentioned analytical approach of Davis and Lauga described the effective slip length 
in dependence on the protrusion angle ip on a surface with grooves perpendicular to 



the flow direction [26[ . The theory assumes rigid bubbles with a full-slip surface, which 
corresponds an infinitely thin liquid-gas interface and vanishing gas density. In this case 
the slip length is given by 

^ = vr(^)/^°°A(.)d., (11) 

with 

s sin 2(f cosh sir + sinh s{tt — 2(f) 



Ais) 



sinh 2s (tt — (/?) + s sin 2(f 



cos 2(p + 



12) 



sinh sir 

The only parameters entering the calculation of the effective slip length b are the 
protrusion angle (p and the ratio between the width of a groove c and the length of 
a unit cell L. The actual values of the slip length have to be calculated numerically. 

In Fig. m we show the dimensionless slip length b/L versus the protrusion angle cp 
as given by Eq. [11] and by the simulation. The system length in this case is L = 40 
lattice units and the width of a groove is c = 30 lattice units. The protrusion angle 
is measured by fitting a circle to the cross section of the interface. Since the interface 
between the vapor and the liquid is diffuse, the actual bubble position is not strictly 
defined. This is in contrast to the theoretical solution and renders the determination 
of the protrusion angle difficult. Therefore, we apply different threshold values gt for 
the fluid density at the interface to determine the protrusion angle. This value has to 
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Figure 4. Slip length versus protrusion angle if for different threshold values. The groove has a 
width of c = 30 and the system length is L — 40. Since the multiphase lattice Boltzmann model 
used is a diffuse interface model, the protrusion angle is not defined uniquely but depends on 
where the actual interface position is assumed, i.e. which threshold gt is chosen. Due to the non 
perfect slip on the vapor phase the analytical solution by Lauga et al. is scaled down by a factor 
of 0.75. 



be somewhere between the high (hquid) density gig = 2.2 and the lower (gas) density 
Qg = 0.05. In Fig. mthe effect of choosing different threshold values between gt = 0.02g\q 
and 0.9piq is demonstrated. Here, it can be observed that a variation of the threshold 
can lead to a shift in the protrusion angle of more than 20°. 

By comparing the simulation data to the results of Davis and Lauga Eq. [HI we 
find that the qualitative shape of the curve is well reproduced, i.e., the slip length first 
increases with rising yj up to a maximum at v^max and then follows a steep decrease for 
high ip. As will be shown later, it eventually can even become negative. However, in 
addition to the possible variation of if we find a second deviation between theory and 
simulation: the detected shp length b is lower than predicted by Davis and Lauga. To 
be able to compare theory and simulation, the theoretical values are scaled by a factor 
of 0.75 to fit the data. This can be explained by the fact that the diffuse interface must 
not be described by a smooth full slip cap. Instead it shows a finite slip due to the 
friction within the interface region. Further, the density ratio between liquid and vapor 
is limited in the lattice Boltzmann model used. In the presented case it is only 1/44. 
Therefore, the shear resistance on the bubble surface is only reduced by this factor, 
while in a real system consisting of, e.g., water and air, this ratio would be of the order 
of 1/1000 rendering the assumption of no shear more realistic. Apart from the shift due 
to the finite interface width, the simulation is able to recover the main conclusions from 
the theory. 

We further investigate the influence of the area covered by the bubble = c/L on 
the slip length. Different channel widths are been investigated for different protrusion 
angles if. The threshold value determining the interface position was chosen to be 
gt = O.OAgiq. Our results are shown in Fig. [5l The qualitative behavior of the results 
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Figure 5. Normalized slip length b/L versus the protrusion angle tp, for different ratios between 
the width of the grooves c and the system length L. The behavior is similar to the prediction of 
Eq. [11] Further, a larger (j> = c/L leads to a strong increase of the slip length. 



again follows Eq. [TTl The maximum slip length is increased very strongly by increasing 
the relative width of the groove. An increase of the surface coverage from = 0.5 to 0.8 
leads to an increase of the dimensionless slip length where the values for the respective 
maxima range from b/L = 0.07 to 0.19. Such a strong increase is also predicted by 
Eq. [m Further, due to the influence of the increased roughness for ip > 60° the slip 
length becomes smaller for larger ip. Also for ip < —40° the slip becomes smaller than 
b < 0.05L and nearly independent on the surface coverage (p. 

To conclude, in the current section we have introduced our simulation methodology 
and demonstrated that it is well suited to reproduce experimental flow profiles and to 
study the influence of the protrusion angle yj on a detected effective slip. However, it 
has to be taken into account that the simulation model is a diffuse interface method 
effectively rendering our bubbles to be on the nanoscale. A comparison with the 
theoretical description of Davis and Lauga demonstrates that the diffuse interface can 
have a strong effect on the actual slip which has to be considered when designing surfaces 
with nanoscale patterns for specific applications. 



3.2. Spherical nanobubbles 

Steinberger et al. utilized surfaces patterned with a square array of cylindrical holes 
to demonstrate that gas bubbles present in the holes may cause a reduced slip [24 . 
Numerically, they found even negative slip lengths for flow over such bubble mattresses, 
i.e., the effective no-slip plane is inside the channel and the bubbles increase the flow 
resistance. Motivated by Steinberger's work, we consider negative slip lengths on bubble 
covered surfaces and also discuss the question of shear-rate dependent slip. In particular, 
we demonstrate that bubbles can generate a shear-rate dependence of the detected slip 
length. The current section is a review of the studies presented in 25|, l58|] and here 
we relate them to recent references in the literature and newer findings as presented in 
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remainder of this contribution. 

The distance between top and bottom walls is d = 40 lattice nodes in all simulations 
in this section and the area fraction of holes is kept at 0.43. Similar to the previous 
section, a unit cell of the regular array is included in a simulation and periodic boundary 
conditions are applied at domain boundaries. 




Figure 6. 3D visualization of a spherical bubble. While the left panel shows the undisturbed 
shape, the right panel demonstrates the deformation of the bubble under strong shear. 



We investigate the effect of a modified protrusion angle and different surface 
patterns by using square, rectangular, and rhombic bubble arrays. The cylindrical holes 
have a radius a = 20 lattice units and the area fraction of the holes is equal in all cases. 
The shear rate is such that the Capillary number Ca = fiaG.s/'y = 0.16. Here, Gg and 7 
are the shear rate and surface tension, respectively. A snapshot of a simulated system 
is shown in the left part of Fig. [61 The slip lengths obtained for differently shaped unit 
cells and protrusion angles are shown in Fig. [71 Following the discussion on the diffuse 
interface and the definition of the protrusion angle as given in the previous section, we 
choose the threshold value for the density defining the interface to be at gt = O.S^iq. 
The observed behavior is similar to that reported in 2J], where a square array of holes 
was studied. In particular, we observe that when (p is large enough b becomes negative. 
Moreover, for the case presented in the current paper the slip length is maximized when 
the protrusion angle is close to zero. In practical applications, this would allow the 
highest possible throughput in a microchannel to be obtained. The behavior of the slip 
length can be explained by thinking of an increased surface roughness if the protrusion 
angle is larger or smaller than zero which is in competition with the increased area of 
low friction in the case of large bubbles. Since the area fraction of the bubbles is kept 
constant for all three different unit cells, our results clearly indicate that slip properties 
of the surface can be tailored not only by changing the protrusion angle but also by the 
array geometry. In the presented study, the highest slip lengths are obtained for the 
rhombic unit cell. 

Generally, our findings are in accordance to the results presented in the previous 



section as it has been confirmed in the theoretical paper of Davis and Lauga 26 



Obviously, also here the limitations of the comparison due to the diffuse interface 
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Figure 7. Slip length 6 as a function of the protrusion angle cj). Corresponding results are given 
for a rhombic array (triangles), a rectangular array (diamonds), and a square array (circles). 



and finite density ratio of the liquid and gas pliases still hold. In addition, the 
simulated system is now three dimensional while Davis and Lauga only considered a 
two-dimensional case. 
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Figure 8. The slip length 6 as a function of the capillary number Ca for a square array of 
spherical bubbles with three different protrusion angles, (p = 63°, 68°, and 71° (from uppermost 
to lowermost) is shown. The inset depicts cross sections of liquid-gas interfaces for four capillary 
numbers and (p — 71° [25|. 



In the following, the shear-rate dependence of the slip length is investigated. As the 
shear rate and thus the viscous stresses increase the bubbles are deformed and the flow 
field is modified. An example of such a deformed bubble is shown in Fig. [6] (right). In 
Fig. [HI we show the simulated normalized slip length b/a as a function of the Capillary 
number for three different protrusion angles. The Capillary numbers chosen are at the 
higher end of the experimentally available range. The inset of Fig. [8] shows cross sections 
of the liquid gas interface of a bubble with protrusion angle (p = 71° demonstrating the 
correlation between shear rate and deformation. Our lattice Boltzmann simulations 
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clearly show shear-rate dependent slip, but the behavior is opposite to that found in 
some experiments: in fact, the slip lengths measured by us decrease with increasing 
shear due to a deformation of the bubbles which effectively reduces the roughness of 
the surface. In the experiments, surface force apparatuses are used (see, e.g., Ref. [l9J), 
where a strong increase in the slip is observed after some critical shear rate. This shear- 
rate dependence has been explained, e.g., with formation and growth of bubbles [ul, l73|. 
In our simulations, there is no formation or growth of the bubbles as we only simulate 
a steady case for given bubbles. The experiments on the other hand are dynamic. 
However, our results indicate that the changes in the flow field which occur due to 
the deformation of the bubbles cannot be an explanation for the shear-rate dependence 



observed in some experiments. Our results are consistent with [6J], where it is shown 
that smaller roughness leads to smaller values of a detected slip. In the present case, 
the shear reduces the average height of the bubbles and thus the average scale of the 
roughness decreases as well. However, a possible explanation of the dependence of the 



slip length on the shear rate was recently given by Gao and Feng [271]. They argue 
that the experimentally observed increase of b can be explained by a depinning of the 
three-phase contact line resulting in the extreme case in a continuous gas film on the 
surface. Such a behavior is omitted in our simulations since the interface between liquid 
and gas favours to stay at the edge of a hole due to the different wettabilities inside the 
holes and on the top of the surface which effectively results in a contact-line pinning. 

3.3. Bubble-size distributions 

As already stated in the introduction, in all previous contributions to the literature 
concerning slip flow over a bubble mattress, only mattresses which consist of evenly 
sized bubbles are taken into account. However, such a situation is not necessarily 
realistic since in reality some kind of imperfections are inevitably present. Therefore, 
this section focuses on mattresses with bubbles of different size. For the sake of simplicity 
we restrict ourselves to two different sizes only. In practice, such a situation can be 
achieved by structuring the surface with holes that have different radii. The height 
of the bubbles occurring on such a structuring can be tuned by variation of the bulk 
pressure. For a given bulk pressure, bubbles located in larger holes protrude stronger 
to the flow channel than those in the smaller ones. We show that the bubble-size 
distribution can significantly affect the shp behavior and, e.g., decrease the negative slip 
caused by strongly protruding bubbles. While this reduction may hamper experimental 
observations of negative slip, it is an attractive candidate for designing bubble surfaces 
aiming to maximal flow slippage. 

The geometrical setup used in the present simulations is shown in Fig. [91 As in the 
previous section we simulate a Couette flow between two parallel plates, where one of 
the plates is smooth and another is patterned with an array of cylindrical holes. The 
holes form two nested square arrays such that one of the arrays consists of holes with 
radius ai and the other one of holes with radius 02- A set of four unit cells of the 
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Figure 9. Visualization of the simulation setup. Four unit cells are shown. Each has small 
bubbles (radius 02) at four edges of the cell connected through periodic boundary conditions, and 
a large bubble at the centre of the cell (radius ai). Thus, the bubble mattress consists of two 
nested square arrays of bubbles. In the case shown 02/01 = 0.67. The upper wall is smooth and 
it is moved with constant speed to achieve a shear flow. 

system is shown in Fig. [9l In all simulations, ai is fixed to 30 lattice units and 02 varies 
between 20 and 30 lattice units. The system size is adjusted accordingly such that the 
area fraction of the holes is = 0.42 regardless of the value of 02- Notice that when 
the radii ai and 02 are unequal, the protrusion angles related to the respective holes are 
also unequal for a given value of the bulk pressure. 
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Figure 10. Slip length as a function of the ratio 02/01. The curves are for different values of 
bulk pressure such that the protrusion angle of the larger bubbles are ip = 80°, 76°, 73°, and 70° 
from bottom to top. The shear rate is fixed and corresponds to capillary number Ca — 0.11. 



First we investigate how a variation of the size of the smaller holes affects the 
effective slip. As the radius 02 is reduced, also the protrusion angle of the bubbles 
attached to the smaller holes decreases, i.e., these bubbles become flatter. In previous 
investigations 24J426l| it has been shown for the case of mono-sized protruding bubbles. 
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that when the protrusion angle decreases, the shp length is increased due to the reduced 
roughness of the bubble mattress. On the basis of these findings, one could expect 
that when the radius a2 is reduced, the slip length decreases despite the fact that the 
radius of the larger holes is kept constant. The simulation results, as shown in Figs. [10] 
and [TT], show that this is indeed the case for small changes in the size of smaller holes. 
However, when the size difference between larger and smaller holes is increased such 
that 02/(21 ~ 0.8, a maximum in the effective slip length is observed and thereafter, for 
even smaller values of 02/01, the slip length starts to decrease. This observation can 
be understood as a competition of two different factors affecting the slip phenomenon. 
The first one is related to the above described flattening of the smaller bubbles, which 
dominates at ratio 02/01 close to unity. However, as the size of the smaller holes is 
further reduced, the area fraction of the smaller bubbles also decreases as well as their 
effect on the total effective slip. On the same basis, one might expect that as 02/01 — )■ 0, 
slip approaches the same value as for O2/01 = 1. However, as the former case corresponds 
to a square array of bubbles and the latter to a rhombic one, the actual values of slip 
length may be different as shown in Ref. 25| or Sec. 13. 2[ 
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Figure 11. Slip length as a function of the ratio 02/01. The curves are for different shear rates 
corresponding to Ca ~ 0.036, 0.11, 0.18, and 0.25 from top to bottom. The bulk pressure is fixed 
such that the protrusion angle of the large bubbles is (p = 70°. 



In Fig. [To] we show how the bulk pressure affects the slip length. One can see 
that when the pressure is increased (i.e., the protrusion angle ip is decreased), also the 
effective slip length increases. This is due to the decreasing overall roughness as the 
height of all bubbles in the system is reduced. From Fig. [10] we also observe that the 
maximum of the slip length becomes less pronounced and shifts towards smaller values 
of ratio 02/01 as the bulk pressure increases. This behavior is reasonable as in this case 
the difference in the heights of larger and smaller bubbles is reduced (see Fig. [T2|) . 

Finally we consider the effect of shear rate on the slip length. It was previously 
shown in section 13.21 that increasing shear rates cause the slip length to decrease. This 
dependence is related to the reduced height of the roughness as the shear deforms the 
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bubbles (see also Ref. 6J]). Similar behavior is found also for the mattress with bimodal 
bubble-size distribution (see Fig. [11]). However, we point out that it has recently been 
shown that depinning of the contact line at high shear rates may lead to formation 
of continuous gas film on the surface, and this depinning transition may change the 
qualitative behavior of the shear-rate dependence 271]. In the present work, contact 
lines are pinned on the hole edges at all shear rates. 

To summarize, in this section we simulated laminar flow past bubble arrays 
consisting of bubbles of two different sizes, and studied the flow slippage caused by the 
bubbles. The slip behavior was found to differ from that previously found for surfaces 
made up of evenly sized bubbles. The observed difference is due to the additional larger- 
scale roughness that results from the more complex surface geometry. The presented 
results indicate that even rather small non-uniformity in the mattress may considerably 
increase the slip and thus complicate the observation of negative slip. On the other 
hand, as the negative slip is usually an unwanted property from the application point of 
view, non-uniformity and disorder could be utilized to prevent negative slip and increase 
the throughput in microfluidic devices. 



4. Conclusion 

To conclude, we presented lattice Boltzmann simulations of laminar flow past structured 
surfaces which are covered with nanobubble arrays. After a qualitative comparison 
to experimentally observed flow over cylindrical bubbles we demonstrated that our 
simulations are able to reproduce the general features predicted by the theoretical 
description of Davis and Lauga. However, for a quantitative comparison the theory 
needs to be extended to take into account diffuse liquid-gas interfaces as well as a flnite 
density difference between liquid and gas. By studying arrays consisting of spherical 
bubbles we demonstrated that a strong shear can deform the bubbles and thus reduce 
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the apparent roughness of the surface. As a consequence the effective shp is reduced. 
We have further shown that patterned surfaces can be designed for maximal shppage 
by variation of the position of the bubbles as well as by utilizing a distribution of their 
size. Our results can be utilized to develop micropatterned superhydrophobic surfaces 
for specific microfluidic applications where a local variation or an especially high value 
of the effective slip might be used to direct flow or to maximize the throughput in a 
microchannel. Further, we believe that our results can contribute to the understanding 
of some controversially discussed experimental results on the shear rate dependence of 
boundary slippage. 
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